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I.  INTRODUCTION 


The  general  field  of  system  modeling  and  parameter  estimation  is  a  rich  and 
current  area  of  research.  In  this  thesis  we  apply  four  principal  concepts  to  solving  the 
least  squares  normal  equations,  which  have  the  general  form: 


ila  =  r  (1.1) 

where  R  is  the  data  correlation  matrix,  a  is  the  parameter  vector  to  be  determined, 
and  r  is  the  cross-correlation  vector.  First,  we  desire  to  iteratively  solve  the  problem 
rather  than  use  direct  inversion  of  the  correlation  matrix  R.  Second,  we  use  a  Toeplitz 
approximation  matrix  splitting  technique  [Ref.  1,  2]  to  set  up  the  iterative  equations. 
Third,  we  increase  the  diagonal  dominance  of  the  correlation  matrix  R  in  order  to 
improve  the  convergence  properties  of  the  iterative  equations.  Finally,  we  partition 
the  norm^d  equations  (1.1)  in  order  to  improve  the  computational  efficiency  and  rate 
of  convergence  of  the  iterative  algorithms.  In  all  these  ceises  our  interest  is  to  study 
both  fixed  data  or  off-line  adgorithms  and  data-adaptive  or  online  algorithms. 

A.  PERFORMANCE  CRITERIA 

The  performance  criterion  for  the  algorithms  developed  in  this  thesis  is  based 
on  a  le£ist  squares  formulation.  The  algorithms  minimize  the  sum  of  the  squared 
error  between  the  true  system  output  and  the  estimated  system  output.  We  choose 
to  use  the  covairizuice  method  for  setting  up  the  problem  in  order  to  remain  within  the 
available  data.  Although  the  covau'iance  method  causes  the  data  correlation  matrix 
R  of  (1.1)  to  be  non- Toeplitz,  it  has  the  potentied  to  provide  a  less  biased  pajameter 
estimate  than  the  autocorrelation  method  [Ref.  3]. 
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B.  THESIS  OVERVIEW 

This  thesis  is  divided  into  five  chapters,  including  the  Introduction.  In  Chapter 
II  we  explicitly  develop  the  least  squares  data  formulation  for  both  the  finite  impulse 
response  (FIR)  and  infinite  impulse  response  (IIR)  system  models.  For  the  HR  system 
we  take  special  care  to  account  for  the  dependency  of  the  filter  output  on  the  input 
and  output  filter  coefficients.  Chapter  III  presents  iterative  algorithms  for  FIR  sys¬ 
tems.  First  we  consider  the  fixed  data  Ccise,  and  present  the  simplistic  Gauss-Siedel 
algorithm  as  a  way  of  introducing  the  concepts  of  matrix  splitting  and  iteration.  Then 
we  present  the  Toeplitz  approximation  matrix  splitti^^g  algorithm,  and  develop  two 
modified  versions  of  it.  Here  we  apply  two  new  techniques:  diagonal  pertubation  and 
multiple  partitioning.  Next,  we  develop  and  contraist  a  new  data-adaptive  algorithm 
with  the  well  known  recursive  least  squares  (RLS)  algorithm.  IIR  systems  are  inves¬ 
tigated  in  Chapter  IV.  First  we  present  the  fixed  data  IIR  Toeplitz  approximation 
iterative  algorithm,  and  then  attempt  to  apply  diagonal  pertubation  and  partition¬ 
ing  to  it.  Achieving  marginal  success  for  the  fixed  data  case,  we  then  turn  to  the 
data-adaptive  case.  Here  we  develop  a  new  data-adaptive  IIR  algorithm  based  on  the 
Toeplitz  approximation  matrix  splitting  that  was  used  in  the  fixed  data  case.  This 
tdgorithm  performs  well,  and  does  not  have  the  stability  problems  eissociated  with 
the  IIR  RLS  algorithm,  which  we  present  for  comparison.  In  the  final  chapter,  we 
summarize  the  results  of  simulation  and  recommend  topics  for  future  research. 
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II.  LEAST  SQUARES  DATA  FORMULATION 

The  central  problem  of  this  thesis  is  to  determine  the  parameters  of  an  optimal 
linear  filter  based  on  the  input  and  output  measurements  of  a  system.  For  the  system 
shown  in  Figure  2.1,  both  the  input  data  sequence  x(n)  and  the  output  data  sequence 
y{n)  axe  known.  The  objective  is  to  determine  a  linear  filter,  or  model  system  which 
produces  the  output  from  the  given  input.  Using  least  squares  minimization  tech¬ 
niques  gives  a  filter  that  is  optimal  when  it  minimizes  the  sum  of  the  squared  error 
between  the  known  output  y{n)  and  the  filter  output  y{n).  The  first  section  of  this 
chapter  contains  the  leaat  squares  data  formulation  for  the  finite  impulse  response 
(FIR),  or  non-recursive  filter.  The  next  section  presents  the  least  squares  data  for¬ 
mulation  for  the  infinite  impulse  response  (HR),  or  recursive  filter,  which  has  some 
subtle  but  significant  differences  from  the  FIR  case.  This  chapter  concludes  with 
some  generalizations  about  the  algorithms  to  be  developed. 


Model 

System 


y(n) 


x(n) 


Unknown 

System 


’  Error 


y(n) 


Figure  2.1:  Fundamental  System  Model 
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A.  FINITE  IMPULSE  RESPONSE  (FIR)  MODELS 

Consider  an  M‘'*-order  FIR  filter  with  input  i(n)  and  output  y{n).  The  output 
at  time  n  is  simply  a  weighted  linear  combination  of  the  input: 

y{n)  =  x{n)a^  +  x{n  —  d - +  x{n  —  M)a^.  (2.1) 

This  can  be  written  as 

y(n)  =  x„a^,  (2.2) 

where x„  =  [a:(n)  z(n  — 1)  •  •  •  x{n  —  M)]  is  a  1  xM+1  input  data  vector,  =  [a^ 

■  ■  ■  is  a  A/  +  1  X  1  vector  containing  the  filter  weights,  or  coefficients,  and  y(n)  is 
the  filter  output  at  time  n.  The  superscript  M  indicates  that  these  are  the  coefficients 
of  an  M*^-order  filter.  Because  the  true  value  of  a^  is  not  known,  equation  (2.2) 
produces  an  estimate  for  the  output,  y{n).  Consequently,  the  error  e(n)  between  the 
true  output  y{n)  and  the  FIR  filter  output  y(n),  is  given  by 

e(n)  =  y(n)  -  y(n)  =  y{n)  -  x„a^.  (2.3) 


The  key  to  a  least  squares  solution  of  equation  (2.2)  for  is  forming  an 


overdetermined  set  of  P  +  1  equations  (where  P  >  M  is  necessaury  for  a  unique 


solution): 


y(«) 

i(n)  x(n  —  1) 

x(n  —  M) 

\a^] 

y(n  -  1) 

x(n  —  1)  x(n  —  2) 

x{n  —  M  —  \) 

af 

.  -  P) . 

x{n  —  P)  x{n  —  P  —  l)  ••• 

x{n  —  M  —  P) 

.  OaJ  . 

Xn 


y 


Xn-l 


(2.4) 


L  Xn-P  J 


this  can  be  written  more  compactly  in  matrix  notation  as 


y  = 


Xa 


M 


(2.5) 
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Equation  (2.5)  is  overdetermined  in  the  sense  that  there  are  more  equations  than 
there  are  unknowns,  or,  equivalently,  there  are  more  rows  of  the  data  matrix  X  than 
there  are  coefficients  in  the  filter  coefficient  vector  a^.  Using  (2.5),  the  error  vector 
can  be  written  as:  e  =  y  —  y  =  y  —  XsJ^.  An  optimal  least  squares  solution  minimizes 
the  sum  of  the  squared  errors  which  is  represented  as: 


c= 

j=n-P 


(y  -  Xsi^fiy  -  Xb^) 

(y^  -  B^'^X'^){y  -  Xb!^) 

y^y  -  y'^XB^  -  B^^X^y  +  b’^^X^Xb.  (2.6) 


The  minimization  is  accomplished  by  taking  the  derivative  of  the  sum  of  the  squared 
errors  ^  with  respect  to  the  filter  coefficient  vector  [Ref.  4]  and  setting  it  equal 
to  zero: 

=  0  -  2X^y  +  2J?’'A:a"  =  0.  (2.7) 

Rearranging  gives 


=  X'^y 


(2.8) 


as  the  requisite  condition  for  to  be  minimum.  It  is  importajit  to  notice  that  the  input 
data  which  forms  X,  and  the  true  output  data  contained  in  y  are  independent  of  the 
filter  coefficients  contained  in  b^ .  Because  (2.8)  contains  the  estimate  of  the  optimal 
M‘^-order  FIR  filter  coefficients  a^,  solving  (2.8)  determines  the  optimal  FIR  filter. 
Therefore,  equation  (2.8)  defines  the  FIR  filter  problem  to  be  solved  in  Chapter  III. 
A  more  standard  representation  of  (2.8),  referred  to  as  the  normal  equation  is 


Rb^ 


=  r, 


(2.9) 


where  R  =  X^X  is  referred  to  as  the  correlation  matrix,  and  r  =  X^y  is  the  cross¬ 
correlation  of  the  input  data  matrix  X  and  the  output  data  vector  y.  Theoretically 
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if  R  is  full  rank,  the  solution  to  (2.9)  is  given  by  =  R~^t.  This  thesis  presents 
iterative  algorithms  which  are  alternatives  to  the  direct  inversion  of  R.  Directly 
inverting  R  to  solve  (2.9)  ha.s  three  major  shortcomings: 

1.  Direct  inversion  is  computationally  intensive;  on  the  order  of  multiplications 
for  a  M  X  M  correlation  matrix  R. 

2.  If  R  is  nearly  singular,  or  very  poorly  conditioned  it  may  not  be  possible  to 
compute  R~^. 

3.  Unlike  iterative  algorithms,  direct  inversion  cannot  provide  an  incremental  es¬ 
timate  of  the  solution. 

B.  INFINITE  IMPULSE  RESPONSE  (HR)  MODELS 

The  general  causal  linear  IIR  filter  is  more  complicated  than  the  FIR  filter 
presented  above.  This  is  because  the  filter  output  y{n)  is  a  linear  combination  of  the 
input  x(n),  x(n  —  !),••  •,  x{n  —  Af),  as  well  as  the  previous  output  y{n  —  1),  y{n  —  2), 

•  •  •)  y{i^  ~  N),  where  M  and  N  are  the  order  of  the  input  and  output  coefficients, 
respectively.  Starting  with  the  difference  equation  representation  of  the  IIR  filter: 

y(n)  =  yin  -  l)b^  +  y{n  -  2)b^  +  ■  •  ■  +  y{n  -  N)b^ 

+  xin)a^  +  x{n  -  l)a^ -\ - \-  x{n  -  M)aM,  (2.10) 

which  can  be  written  as 

=  [  yn-i  :  x„  ] 

y(n)  =  ZnO,  (2.11) 

where  z„  =  [y„_i  ;  Xn]  =  [y(n  -  1)  y(n  -  2)  •  ■  •  y(n  -  N)  :  x{n)  z(n  -  1)  •  •  ■  x{n  -  M)] 
isal  xiV-fAf-fl  data  vector,  9  =  ^  ^  •' 


6 


isaiV  +  M  +  1  xl  vector  containing  the  HR  filter  weights.  Again,  because  the  true 
value  of  and  are  not  known,  equation  (2.11)  will  only  produce  an  estimate 
for  the  output,  y(n).  Thus,  the  error  between  the  true  output  y{n)  and  the  HR  filter 
output  j/(n),  is  given  by 


e(n)  =  y{n)  -  y{n)  -  y{n)  -  z„e. 


(2.12) 


As  before,  the  least  squares  solution  to  (2.11)  for  6  comes  from  an  overdeter¬ 
mined  set  of  P  1  equations  (where  P  >  N  +  M  is  necessary  for  a  unique  solution); 

y{n) 

y{n-l) 


y{n  -  P)  \ 

y{n  -  1) 
y{n  -  2) 


y{n  —  N)  :  i(n) 

y{n  —  iV  —  1)  :  x{n  —  1) 


x(n  —  M) 
x(n  —  M  —  1) 


y(n  —  P  —  1)  y(n  ~  N  —  P)  :  x{n  —  P)  x{n  —  M  —  P) 


y  = 


Yn-l 


Xn 

Zn 

Xn-1 

e  = 

Zn-1 

Xn-P 

.  2n-P  . 

nM 

J 


(2.13) 


or,  in  matrix  notation 


y  =  ze. 


(2.14) 


where  the  data  matrix  Z  =  [z„  z„_i  ■  •  •  which  is  used  to  form  am  error  vector. 


e  =  y  ~  Z6.  Once  again,  the  optimad  least  squaires  solution  minimizes  the  sum  of  the 
squaured  errors  which  now  has  the  form: 


(2.15) 


C= 

j=n-P 

=  {y  -  Ze)'^{y  -  Ze) 

=  {y^-e'^z'^){y-Ze) 

=  y^y-y'^Ze  -e'^Z^y  +  e^Z'^Ze. 


If  the  gradient  terms  that  arise  from  formal  differentiation  are  neglected,  the  derivative 
of  the  sum  of  the  squared  errors  (  with  respect  to  the  filter  coefficient  vector  0  is: 


m 

8(0) 


0  -  2Z'^y  +  2Z'^Z0  =  0. 


(2.16) 


Rearranging  gives 


z'^ze  =  Z'^y 


(2.17) 


as  the  requisite  condition  for  (  to  be  minimum.  Again,  (2.17)  contains  the  estimate 
of  the  HR  filter  coefficients  0 ,  and  its  solution  determines  the  optimal  least  squares 
HR  filter.  Therefore,  equation  (2.17)  defines  the  HR  filter  problem  to  be  solved  in 
Chapter  IV,  and  the  normal  equation  representation  is 


R0  =  r, 


(2,18) 


where  R  =  Z^Z  is  referred  to  as  the  correlation  matrix,  and  r  =  Z^y  is  the  cross¬ 
correlation  vector. 

The  differentiation  in  (2.16)  appears  quite  simple  because  it  ignores  the  de¬ 
pendence  of  the  filter  output  y{n)  on  the  previous  filter  input  and  output.  Formal 
differentiation  of  ^  gives  [Ref.  5]: 


m 

d{0) 


T,j=n~piyij)  -  yU)T 

.  alw  Ej=n-p(i/(i)  -  y{j)Y 


=  0. 


(2.19) 
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Letting  and  q:„  =  represent  the  vector  partial  derivatives  of  tne  filter 

output,  we  can  write 

E  f  f'  1(!'(j)-M)  =  0,  (2.20) 

j=n-P 

or,  rearranging  and  letting  rl^j  =  [^J  gives 

V’jZj 

p 

(2.21) 

R  ^ 

so  that  we  have  the  same  form  as  the  normal  equations  of  (2.18),  R6  =  r.  However, 
R  no  longer  represents  the  autocorrelation  of  the  data  matrix  Z.  In  this  case  R 
represents  the  cross- correlation  matrix  between  the  gradient  components  of  the  filter 
output  ^  and  the  filter  data  matrix  Z,  and  r  represents  the  cross-correlation  between 
^  and  the  true  output  vector  y. 

Regardless  of  the  formulation,  a  standard  solution  to  (2.18)  is  given  by  0  = 
R~^v.  Again,  the  objective  is  to  present  iterative  algorithms  which  are  alternatives 
to  the  direct  inversion  of  R.  Finally,  for  the  algorithms  presented  in  this  thesis,  there 
are  three  general  characteristics  governing  their  suitability: 

1.  The  algorithm  should  use  the  minimum  amount  of  data  required  to  produce  a 
solution. 

2.  The  £ilgorithm  should  converge  quickly.  The  desired  number  of  iterations  of  the 
adgorithm  should  not  exceed  2  to  3  times  the  number  of  parameters  being  solved 
for. 

3.  The  computational  complexity  of  the  algorithm  should  be  less  th2ui  that  of 


direct  inversion. 


III.  FINITE  IMPULSE  RESPONSE  (FIR) 

SYSTEMS 

In  this  chapter  we  consider  the  finite  impulse  response  (FIR)  system.  We  develop 
two  general  categories  of  solution  algorithms,  fixed  data  and  data-adaptive.  The  fixed 
data  algorithms  are  the  Gauss-Seidel  iterative  method  and  the  Toeplitz  approxima¬ 
tion  iterative  method.  The  data-adaptive  algorithms  are  the  recursive  least  squares 
(RLS)  method  and  Toeplitz  approximation  method.  Before  developing  the  solution 
algorithms,  it  is  worth  noting  that  the  basic  problem  formulation  described  in  Chap¬ 
ter  II  is  a  covariance  rather  than  a  correlation  formulation  of  the  data;  therefore,  the 
matrix  R  in  =  r  is  symmetric  but  not  Toeplitz.  Accordingly,  we  cannot  apply 
the  Levinson  recursion  algorithm  [Ref.  6]  which  requires  Toeplitz  structure  in  the 
correlation  matrix.  In  our  formulation,  the  data  matrix  is  formed  using  the  covari¬ 
ance  method,  which  has  the  potential  to  provide  a  less  biased  least  squares  solution 
than  a  correlation  method  formulation  [Ref.  3]. 

A.  FIXED  DATA  ALGORITHMS 
1.  Gauss-Seidel  Method 

A  very  simple  and  straightforward  iterative  adgorithm  is  the  Gauss-Seidel 
method  [Ref.  7].  We  drop  the  superscript  M  from  for  simplicity.  Unless  otherwise 
stated,  we  are  considering  an  M‘'‘-order  FIR  system.  Starting  with  (2.9),  Ra  =  r, 
split  R  into  L  +  D  U,  where  L  is  a  matrix  containing  the  strictly  lower  triangular 
elements  of  R,  C/  is  a  matrix  containing  the  strictly  upper  triangular  elements  of  R, 
and  D  is  a  diagonal  matrix  containing  the  main  diagonal  elements  of  R.  Substituting 


10 


this  into  (2.9)  gives 


(I  +  D  +  C/)a  =  r,  (3.1) 

or 

(X  +  X))a  =  -Usi  +  r.  (3.2) 

Making  (3.2)  into  an  iterative  algorithm  requires  that  we  isolate  the  parameter  vector 
a  on  the  left  side  of  the  equation.  Pre-multiplying  both  sides  of  the  equation  by 
(X  +  X»)-^: 

(X  +  X>)-'(X  +  X))a  =  (X  +  D)-^{-Ua  +  r)  (3.3) 

we  have 

a  = -(X  +  D)-'t/a  +  (X  +  X))-'r.  (3.4) 

Note  that  (L  +  D)  must  be  nonsinglaj  for  this  technique  to  work.  Now  there  are  two 
occurrences  of  the  parameter  vector  a.  We  define  the  one  on  the  right  side  of  (3.4)  as 
the  current  value  a^*^  and  the  one  on  the  left  side  as  the  updated  value  where 

k  is  the  iteration  index.  Thus,  the  final  form  of  the  Gauss-Seidel  iterative  aJgorithm 
is: 

=  -(X  +  X>)-^C/aW  +  (X  +  X>)*'r.  (3.5) 

Unfortunately,  this  algorithm  takes  an  excessive  number  of  iterations  to  converge  to 
the  true  parameter  vedues  when  R  is  even  moderately  ill-conditioned.  In  an  effort 
to  improve  this  performauice,  we  applied  the  successive  overrelajcation  (SOR)  accel¬ 
eration  technique  [Ref.  9]  to  (3.5).  The  SOR  technique  requires  that  we  add  an 
acceleration  parameter  u;  as  follows: 

=  -(u;X  -t-  D)-^  [(1  -  u;)Z?  -  +  {uL  -f-  Z?)-'u;r.  (3.6) 

Notice  that  when  u)  equals  1,  (3.6)  degenerates  into  the  stand2U’d  Gauss-Seidel  algo¬ 
rithm.  In  order  to  speed  up  the  rate  of  convergence  of  the  algorithm,  the  acceleration 
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parameter  w  must  assume  its  optimum  value,  which  is  given  by  [Ref.  9]: 


<^opt 


2(1  -  yi  -  /'»«) 


r^max 


(3.7) 


In  (3.7)  the  term  fimax  represents  the  largest  eigenvalue  of  the  Jacobi  matrix,  which 
has  the  form  D~^{—L  —  U).  Unfortunately,  though  this  is  a  popular  and  usually 
successful  acceleration  technique,  for  the  basic  problem  of  this  thesis  it  was  found  to 
produce  little  improvement  in  the  performance  of  the  Gauss-Seidel  algorithm. 

In  order  to  observe  the  performance  of  the  algorithms  developed  in  this 
chapter,  they  are  simulated  with  the  following  FIR  system: 


y{n)  =  0.5x(n-l)  +  0.25a:(n-2)-0.5x(n-4)  +  0.053x(n-5)-|-  0.0345x(n-6) 
-  0.76x(n  -  7)  -f  3.5x(n  -  8)  -  1.0032x(n  -  9)  -  0.0031x(n  -  10).  (3.8) 

Figure  3.1  shows  the  generally  slow  and  gradual  rate  of  convergence  of  the  Gauss- 
Siedel  and  Gauss-Siedel  with  SOR  algorithms.  Each  parameter  track  line  of  Figure 
3.1  corresponds  to  a  coefficient  in  (3.8).  Notice  that  as  the  number  of  iterations  k 
increcises,  the  parameter  track  lines  flatten  out  and  approach  the  values  in  (3.8). 

2.  Toeplitz  Approximation  Algorithm 

One  reason  for  the  slow  rate  of  convergence  of  the  Gauss-Seidel  method  is 
that  it  fails  to  capitalize  on  the  structure  of  the  covariance  matrix  R.  The  covariance 
matrix  is  symmetric  but  not  Toeplitz.  However,  it  does  have  an  underlying  Toeplitz 
structure.  This  point  becomes  clear  when  the  number  of  data  points  used  to  form 
the  matrix  is  increjised;  the  covariance  matrix  then  approaches  a  Toeplitz  structure. 
Toeplitz  matricies  have  many  nice  properties,  and  efficient  methods  such  as  the  Levin¬ 
son  recursion  are  available  to  invert  a  Toeplitz  matrix.  An  iterative  algorithm  which 
takes  advantage  of  the  near  Toeplitz  structure  of  the  covariance  matrix  [Ref.  1]  is  de¬ 
veloped  by  splitting  R  into  a  Toeplitz  matrix  T,  and  a  residual  matrix  5.  The  matrix 
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Figure  3.1:  Parameter  tracks  and  general  rate  of  convergence  for  the  stan¬ 
dard  FIR  Gauss-Seidel  and  Gauss— Seidel  with  SOR  algorithms. 

T  is  obtained  by  averaging  the  diagonal  elements  of  R.  This  Toeplitz  approximation 
provides  a  natural  splitting  of  the  covariance  matrix,  R  =  T  +  S,  which  is  used  to 
develop  an  iterative  algorithm  as  follows.  Beginning  with 


Ra  =  r, 


substituting  for  R  =  T  -I-  5  gives 

(T-h5)a  =  r, 


(3.9) 


(3.10) 


or,  rearranging  we  have. 


Ta  =  r  —  5a. 


But  S  =  R  —  T,  thus 


Ta  =  T  —  {R  —  r)a 

and  isolating  a  gives 

a  =  r-^r-r-^i2a-ha. 


(3.11) 


(3.12) 


(3.13) 
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Finally,  the  iterative  algorithm  has  the  form: 

a(*+i)  =  (3.14) 

And  if  we  let  ao  =  T“^r  be  the  initial  estimate  of  the  parameter  vector  a,  then  we 
have: 

=  ao  +  a<^)  -  (3.15) 

The  Toeplitz  approximation  matrix  T  can  also  be  viewed  in  another  light 
[Ref.  8].  It  can  be  thought  of  as  a  pre-conditioning  matrix  for  R.  Indeed,  part 
of  the  difficulty  with  the  Gauss-Seidel  method  of  the  previous  section  is  that  as  R 
becomes  more  ill-conditioned,  the  iterative  algorithm  takes  i,  uch  longer  to  converge. 
We  observed  that  the  product  T~^R  does  in  fact  have  a  better  condition  number  than 
R  by  at  least  an  order  of  magnitude  when  R  is  very  ill-conditioned.  Thus,  we  might 
expect  that  the  algorithm  of  (3.15)  will  converge  to  the  true  parameter  vector  a  faster 
than  (3.5)  or  (3.6).  This  is  in  fact  the  case.  In  Figure  3.2  it  is  clear  that  the  Toeplitz 
approximation  algorithm  gives  much  faster  convergence  to  the  true  parameter  values 
than  the  Gauss-Seidel  algorithm. 

3.  Toeplitz  Approximation  with  Diagonal  Pertubation  Algorithm 
One  significant  shortcoming  of  (3.15)  is  that  it  fails  to  converge  as  the 
covariance  matrix  R  becomes  extremely  ill-conditioned.  Despite  the  fact  that  T 
functions  as  a  pre-conditioner  for  R,  the  algorithm  of  (3.15)  will  not  converge  as  the 
spectral  radius  (or  largest  eigenvaJue)  of  {T~^S)  becomes  greater  than  one  [Ref.  9]. 
Generally,  the  condition  number  of  R  is  an  indication  of  now  poorly  conditioned  it 
is.  We  observe  that  ^ls  the  number  of  data  points  used  to  form  R  is  reduced  from 
a  large  number  (i.e.,  5  to  10  times  the  order  of  the  filter  being  solved  for)  to  near 
the  minimum  number  required  for  a  unique  solution  (i.e.,  the  order  of  the  filter),  the 
condition  number  of  R  grows  quite  large.  Under  these  circumstances  the  algorithm  of 
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Gauss-Siedel.  P=50 


Toeo  Anrox.  P=50 


Figure  3.2:  Parameter  tracks  and  general  rate  of  convergence  for  the  FIR 
Toeplitz  approximation  and  the  Gauss-Seidel  algorithms. 

(3.15)  does  not  converge.  By  constrast  the  Gauss-Seidel  algorithm  always  converges 
to  the  true  parameter  values,  although  it  takes  an  excessivly  large  number  of  iterations 
for  the  parameter  values  to  approach  their  true  values  (see  Figure  3.3).  Viewing  the 
ill-conditioned  nature  of  i?  as  an  obstacle  to  satisfactory  performance  of  (3.15),  we 
now  introduce  a  technique  which  improves  the  condition  of  R  ajid  permits  a  modified 
version  of  (3.15)  to  converge  as  the  number  of  data  points  approaches  the  minimum 
required. 

Beginning  with  the  fact  that  any  identity  matrix  I  is  invertible  with  a 
condition  number  of  one,  we  seek  to  alter  R  such  that  it  becomes  more  like  aji  identity 
matrix.  This  is  accomplished  by  adding  a  diagonal  matrix  to  R  such  that  the  resulting 
matrix  has  its  main  diagonal  elements  at  least  am  order  of  magnitude  larger  than  any 
other  element  in  the  matrix,  then  this  new  matrix  has  increased  diagonad  dominance 
aind  begins  to  look  more  like  an  identity  matrix.  The  consequence  of  applying  this 
technique  is  that  we  improve  the  condition  of  R  and  the  performance  of  (3.15). 


Figure  3.3:  Parameter  tracks  and  general  rate  of  convergence  for  the  FIR 
Toeplitz  approximation  and  Gauss-Seidel  algorithms. 

If  we  let  (T  be  a  scalar  multiplier  and  the  Toeplitz  approximation  splitting 
of  R  be  defined  as  R  =  To  -f  5,  then  we  can  derive  a  modified  version  of  (3.15)  as 
follows.  Beginning  with 


=  r,  (3.16) 

adding  a  diagonal  matrix  (scaled  by  <7)  to  R  gives 

(R  +  cr/)a  =  r  +  era  (3. 17) 

and  substituting  for  R  =  To  +  S  we  have 

(To  +  5  +  (7/)a  =  r  +  era.  (3.18) 

Now,  letting  T  =  To  +  <tI,  we  can  write 

Ta  =  r  +  <Ta-5a  (3.19) 

but  S  —  R  +  <t1—  T,  and  substituting  gives 

Ta  =  r  +  <Ta- (R  +  er/- T)a.  (3.20) 
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After  rearranging  it  is  easy  to  see  that  the  aa.  terms  cancel 

Ta  =  r  +  <Ta  —  era  —  (i?  —  T)a  (3.21) 

so  that  we  can  now  write 

a  =  r-^r-r-'i?a  +  a.  (3.22) 

Now,  letting  ao  =  T~^r  and  rearranging,  the  final  form  of  the  Toeplitz  approximation 
with  diagonal  pertubation  iterative  algorithm  is: 

=  ao  +  .  (3.23) 

Interestingly,  the  algorithm  of  (3.23)  has  exactly  the  same  form  as  (3.15).  The  only 
difference  is  that  the  matrix  T  in  (3.23)  is  the  original  Toeplitz  approximation  matrix 
To  plus  some  diagonal  matrix  <jI.  The  new  algorithm  of  (3.23)  was  observed  to 
converge  where  the  algorithm  of  (3.15)  fails.  It  always  reaches  the  vicinity  of  the  true 
parameter  values  quickly,  but  requires  a  lairge  number  of  iterations  to  converge  to 
the  exact  solution  (see  Figure  3.4).  In  simulation,  we  observe  that  optimum  vcilues 
for  <T  axe  inversely  proportionaJ  to  the  condition  number  of  R.  One  difficulty  with 
this  algorithm  is  that  determining  the  optimum  value  for  cr  is  a  problem  as  it  is  data 
dependent.  Therefore,  achieving  optimal  performance  of  this  algorithm  is  dependent 
on  finding  a  suitable  value  for  <r.  We  were  unable  to  find  a  general  expression  for  the 
optimum  value  of  a  in  this  thesis. 

4.  Partitioned  Toeplitz  Approximation  Algorithm 

An  aJtemative  to  improving  the  condition  of  R  in  (3.15)  is  to  block  partition 
R.  It  may  be  partitioned  into  equally  or  unequally  sized  parts,  however  the  number 
of  row  partitions  must  equaJ  the  number  of  column  p^u:titions.  We  now  derive  a 
partitioned  version  of  (3.15).  Beginning  with 

Ra  =  r,  (3.24) 
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Figure  3.4:  Parameter  tracks  and  general  rate  of  convergence  for  the  FIR 
Toeplitz  approximation  and  Toeplitz  approximation  with  diagonal  pertu- 
bation  algorithms. 


partitioning  R  into  four  parts  gives 


Rn 

:  Ri2 

. 

ai' 

ri' 

R21 

.^*2. 

and  substituting  for  R\\  =  Tn  +  Su  and  R22  =  T22  +  S22  we  have 


Til  +  5ii 

•0 

_ 1 

31' 

ri' 

R2\ 

•  T22  +  522  , 

.^2. 

.**2. 

Multiplying  gives  the  following  coupled  equations: 

(Til  +  511)81  +  R12&2  =  Ti 

(T22  +  522)32  +  R2\^\  —  r2 


(3.25) 


(3.26) 


(3.27) 


and  rearranging  this  yields 


Tiiai  =  ri  —  R\2^2  ~  5iiai 
T2232  ~  —  -^2131  —  52232- 


(3.28) 


But  Sn  =  Rii  —  Til  and  S22  =  R22  —  T22,  thus  we  have 


ai  —  T^i^Ti  —  T^i^{Ri2&2  +  (-R11  ~  Tii)&i) 

^2  =  T^^r2  —  T2i^(R2iAi  +  (R22  —  222)^2)  (3.29) 

and  as  before  we  let  aio  =  Ti'i^  ri  and  a2o  =  T22'^2  so  that  now  we  can  write 

ai  =  aio  +  Ui  —  rjj^(i?i2a2  +  2?iiai) 

^2  =  ^20  +  ^2  —  T^^(R21&i  +  R22^2)-  (3.30) 

The  final  form  of  the  partitioned  Toeplitz  approximation  iterative  algorithm  is 

.  =  aio  +  -  rn'(i?i2ai'')  +  i?naS"^) 

=  a2o  +  a(">-r22'(i22iaS"'''^  +  i?22ai"^  (3.31) 

This  partitioning  algorithm  was  also  observed  to  converge  to  the  true  solu¬ 
tion  where  the  algorithm  of  (3.15)  fails.  Additionally,  this  algorithm  converges  to  the 
true  solution  in  about  the  same  number  of  iterations  as  that  of  (3.23),  and  this  algo¬ 
rithm  does  not  require  the  choice  of  an  optimal  parameter  value  such  as  cr  in  (3.23). 
It  is  possible  to  derive  iterative  algorithms  using  this  scime  methodology  beyond  the 
two  pau-titions  shown  in  (3.25),  however  we  did  not  observe  a  significant  difference  in 
performance  of  these  algorithms  over  the  one  shown  in  (3.31)  (see  Figure  3.5).  Fi¬ 
nally,  it  is  also  possible  to  combine  diagonal  pertubation  and  partitioning  in  a  single 
algorithm,  but  again  we  did  not  observe  any  signific2uit  improvement  in  performance. 


B.  DATA-ADAPTIVE  ALGORITHMS 

Although  the  fixed  data  2dgorithms  above  do  provide  a  solution  to  the  least 
squares  problem  i?a  =  r,  they  require  that  a  minimum  amount  of  data  be  accumulated 
before  the  iterative  zdgorithm  may  begin  and  a  solution  determined.  Consequently, 
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Figure  3.5:  Parameter  tracks  and  general  rate  of  convergence  for  the  FIR 
partitioned  Toeplitz  approximation  and  double-partitioned  Toeplitz  ap¬ 
proximation  algorithms. 

real-time  processing  of  signals  as  they  axe  available  is  not  possible.  Data-adaptive 
algorithms  compute  a  solution  based  upon  the  inclusion  of  new  data  with  each  iter¬ 
ation.  They  provide  a  means  of  providing  a  real-time  system  identification  solution. 
In  this  section  we  develop  a  data-adaptive  algorithm  by  extending  the  Toeplitz  ap¬ 
proximation  technique  presented  in  the  previous  section.  Also,  we  present  the  well 
known  RLS  algorithm  for  comparison. 

1.  Toeplitz  Approximation  Algorithm 

In  order  to  avoid  confusion  with  previous  developments  of  the  basic  problem 
for  the  fixed  data  case,  we  define  the  following: 


x„  =  [x(n)  x(n  —  1)  •••  x(n  — M)]^ 

■  oo  ■ 

Ol 

a  = 

-Om- 


(3.32) 

(3.33) 
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so  that  now  y(n)  =  for  a  specific  value  of  the  output  of  aji  M‘^-order  filter.  For 

the  generad  problem  development,  we  have  the  following; 

x(n)  x(n  —  1)  •••  x{n  —  M) 

i(n  —  1)  x(n  —  2)  •••  i(n  — M  — 1) 


i(n  —  P)  x(n  —  P  —  1) 
or,  in  vector  notation  we  have; 


x(n  —  P  —  M) 


Oq 

y(^) 

Oi 

y(n  -  1) 

.OAf . 

,y(n-P) 

(3.34) 


‘•n-1 


cl- 


a  =  y, 


(3.35) 


and  now  it  is  easy  to  see  that  the  baisic  problem  has  the  form; 


,  Xn  Xn— ll  3Cn— p 


"n-l 


lK-pj 


a  =  [xn  x„_i  •••  x„_p]y. 

' - ^ - ' 

r 


(3.36) 


R 


It  is  important  to  realize  how  the  data-adaptive  algorithm  can  be  derived  from  R  and 
r.  Using  a  subscript  n,  we  can  express  R  and  r  in  time  varying  notation  as; 

Rn  =  12 

j=n—P 

=  Xn-pX^_p  +  X„_p+iX^_p+l  +  •  •  •  +  Xn_iX^_l  +X„X^ 

^  '  '  ^1— ^  ■  II  ,  ✓ 


■^^n— 1 


Rn  =  Rn-1  +  X„x: 


(3.37) 


amd 


Tn  =  X)  ^Mj) 

j=n-P 

=  Xn-py(n  -  P)  +  Xn_p+iy(n  -  P  +  1)  +  •  •  •  +  x„_iy(n  -  1)  +Xny(n) 


r„_i 


Fn  =  r„_i  +  Xny(n).  (3.38) 

Now  that  we  have  Rn  and  r„,  we  cam  form  the  data-adaptive  Toeplitz 
approximation  algorithm  ais  in  the  fixed  data  caise; 


Rfi^  — 


(3.39) 
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substituting  Ji„  =  T„  +  Sn  yields 


(T„  +  Sn)a  =  r„,  (3.40) 

cuid  rearranging  gives 

T„a  =  r„  -  ^na. 

But  Sn  —  Rn  —Tn,  thuS 

Tn^  =  Tn  ~  {Rn  ~  ^n)a 
and  isolating  the  parameter  vector  a  we  have 

a  =  -  T~^Rna  +  a 

Finally,  the  iterative  algorithm  has  the  form: 

a(fe+i)  =  +  r-'r„  -  T'^RnsS’^l  (3.44) 

This  algorithm  converges  quickly  to  the  true  parameter  values  (see  Fig  3.6). 
Notice  that  both  2dgorithms  begin  with  time  n  =  fc  =  0.  The  data-adaptive  adgorithm 
begins  to  estimate  parameter  v<dues  immediately,  but  the  fixed  data  algorithm  must 
wait  until  a  minimum  of  Af + 1  data  points  have  been  accumulated  before  it  can  begin 
estimating  parameter  values  for  an  M**'Order  system.  Although  the  data-adaptive 
algorithm  (3.44)  outperforms  the  fixed  data  <dgorithm  of  (3.31),  there  is  a  considerable 
computational  tradeoff.  The  matricies  T~^  and  Rn,  and  the  vector  r„  must  be  updated 
at  each  iteration  of  the  data-adaptive  algorithm.  The  matrix  T~^  is  updated  using 
the  matrix  inversion  lemma  [Ref.  10],  and  Rn  amd  r„  are  updated  from  (3.37)  and 
(3.38),  respectively.  This  requires  approximately  5M’  additional  multiplications  for 
each  iteration  of  (3.44)  as  compared  with  the  fixed  data  adgorithm. 


(3.41) 

(3.42) 

(3.43) 
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Figure  3.6:  Parameter  tracks  and  general  rate  of  convergence  for  the  FIR 
data-adaptive  and  fixed  data  Toeplitz  approximation  algorithms. 

2.  Recursive  Least  Squares  (RLS)  Algorithm 

In  order  to  evaluate  the  effectiveness  of  the  data-adaptive  Toeplitz  approx¬ 
imation  algorithm  above,  it  must  be  compared  against  some  benchmark.  The  RLS 
algorithm  is  recognized  as  a  particularly  optimail  approach  to  the  least  squares  FIR 
problem.  A  simple  presentation  of  the  update  algorithm  is  [Ref.  5]: 

a(*+i)  =  -h  aR;|iX„c(n)  (3.45) 

where  a  is  referred  to  as  the  forgetting  factor.  A  performance  comparison  between  this 
algorithm  and  that  of  (3.44)  reveals  the  astonishing  fact  that  they  perform  exactly 
the  same.  Apppendix  A  shows  several  examples  of  this  under  different  conditions  of 
noise,  and  over-  and  under-modeling.  On  close  analysis,  we  note  that  the  RLS  algo¬ 
rithm  (3.45)  can  be  expanded  and  rewritten  in  the  saune  form  as  (3.44).  Although 
the  algorithms  look  and  perform  the  same,  they  au^e  in  fact  different.  Notice  that  Rn 
jind  Tn  are  summed  over  the  current  and  past  data  in  (3.44),  whereas  these  terms  in 
the  RLS  algorithm  a^e  composed  of  only  the  current  data  values  {Rn  =  x^x„  and 
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Figure  3.7:  Parameter  tracks  and  general  rate  of  convergence  for  the  FIR 
RLS  and  data-adaptive  Toeplitz  approximation  algorithms. 

r„  =  x„t/(n)).  Interestingly,  as  shown  in  Figure  3.7,  the  performance  of  the  algo¬ 
rithms  is  still  the  same.  Unfortunately,  computing  a  more  complicated  R,  and  Tn  in 
the  Toeplitz  Approximation  data-adaptive  algorithm  requires  an  increase  of  approx¬ 
imately  2M^  —  M  multiplications  over  RLS.  Although  this  increased  computation 
does  not  produce  any  improvement  in  performance  over  FIR  RLS  systems  modeling, 
it  will  give  improved  performance  over  HR  RLS,  as  we  shall  see  in  Chapter  IV. 


24 


IV.  INFINITE  IMPULSE  RESPONSE  (UR) 

SYSTEMS 

In  this  chapter  we  consider  the  infinite  impulse  response  (HR)  system.  Here 
again  we  develop  two  general  categories  of  solution  algorithms,  fixed  data  and  data- 
adaptive.  The  slow  rate  of  convergence  of  the  Gauss-Seidel  iterative  algorithm  applies 
to  HR  systems  in  the  same  way  it  applied  to  FIR  systems  in  Chapter  III.  Therefore, 
in  this  chapter  we  only  present  one  fixed  data  adgorithm,  the  Toeplitz  approximation 
iterative  method.  The  data-adaptive  algorithms  are  the  recursive  least  squares  (RLS) 
method  and  Toeplitz  approximation  method.  The  same  discussion  regarding  the  use 
of  techniques  that  are  Toeplitz  dependent  from  Chapter  III  is  applicable  here  as  well; 
our  covariance  formulation  of  the  problem  precludes  us  from  using  such  techniques. 

A.  FIXED  DATA  ALGORITHMS 

1.  Toeplitz  Approximation  Algorithm 

Here  we  present  a  Toeplitz  approximation  iterative  algorithm  [Ref.  2]  of 
the  scime  structure  as  in  the  FIR  case.  One  major  difference  is  that  we  will  start  with 
a  partitioned  version  of  the  algorithm  similar  to  that  of  (3.31).  This  is  quite  natural 
since  we  are  solving  for  two  parameter  vectors.  One  vector,  a^,  represents  the  filter 
input  or  feedforward  coefficients  and  the  other,  b^,  represents  the  filter  output  or 
feedback  coefficients.  As  in  Chapter  III,  we  drop  the  superscripts  M  and  N.  Unless 
otherwise  stated,  the  input  coefficient  vector  has  order  M,  and  the  output  coefficient 
vector  has  order  N.  Thus,  from  Chapter  II,  (2.18),  we  begin  with 

Re  =T  (4.1) 
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and  partitioning  R  gives 


Ryy 

:  Ryx 

b' 

’ry' 

_Rxy 

:  Rxx . 

a 

.Tr. 

We  now  split  the  matricies  Ryy  and  Rxx  into  a  Toeplitz  and  a  residual  matrix: 


Tyy  +  Syy 

H 

_ 1 

b' 

'ry' 

Rxy 

:  Txx  +  Sxx  . 

a_ 

Multiplication  gives  the  following  coupled  equations: 


(4.2) 


(4.3) 


(^yv  ‘^yy)^  "h  —  '^y 

(Tix  +  5rr)a  +  Rxyh  =  Tx  (4.4) 


and  rearranging  these  gives 


Tyyh  -  Ty-  Ryx^  “  ‘^W  ^ 

Txapa  *  Vx  "  *^xxa.  (4.o) 


But  Syy  =  Ryy  ~  Tyy  BJld  SxX  =  RxX  "  tfaUS 

b  =  Tyy  Vy  —  Tyy  {Ryx&  +  {Ryy  ~  ryy)b) 

a  =  r;;^rx  -  r;;*(jRx„b  +  (i?xx  -  rxx)a)  (4.6) 


amd  letting  bo  =  T~^ry  and  a©  =  T~^Tx  as  the  initiad  estimates  of  the  parameter 
values,  we  have 


b  =  bo  +  b  —  T^{Ryx^  +  Ryy^) 

a  =  ao  +  a-r;;^(i2xyb  +  i2xxa).  (4.7) 

The  finaJ  form  of  the  fixed  data  HR  Toeplitz  approximation  iterative  algorithm  is: 


Iterations,  k 
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Figure  4.1:  Parameter  tracks  and  general  rate  of  convergence  for  the  HR 
Toeplitz  approximation  algorithm. 

=  bo  + 

a(*+')  =  ao  +  aW_T,-'(i?^^b(*+^)  +  /2„a('')).  (4.8) 

In  order  to  evaluate  the  general  performance  of  the  algorithms  developed 
in  this  chapter,  the  following  HR  system  is  used  for  simulation: 

y{n)  =  2.7600j/(n  —  1)  —  3.8090y(n  —  2)  +  2.6540y(n  —  3)  —  0.9240y(n  —  4) 

+  i(n)-0.9x(n- l)T0.81x(n-2)-0.65x(n-3)  +  0.36z(n-4).  (4.9) 

Figure  4.1  shows  that  this  algorithm  gives  performance  similar  to  that  of  the  FIR 
adgorithm  (3.15).  It  converges  relatively  quickly  when  R  is  well-conditioned,  amd  fadls 
to  converge  to  the  true  pareuneter  vaJues  as  R  becomes  ill-conditioned. 

2.  Toeplitz  Approximation  with  Diagonal  Pertubation  Algorithm 
As  in  the  FIR  case  a  significant  shortcoming  of  (4.8)  is  that  it  fails  to 
converge  as  the  covariance  matrix  R  becomes  extremely  ill-conditioned.  We  again 
observe  that  as  the  number  of  data  points  used  to  formulate  R  is  reduced  from  a 
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large  number  (i.e.,  5  to  10  times  the  order  of  the  HR  filter  beinc^  solved  for)  to  near 
the  minimum  number  required  for  a  unique  solution  (i.e.,  the  order  of  the  filter),  the 
condition  number  of  R  grows  quite  large.  Under  theoe  circumstances  the  algorithm 
of  (4.8)  does  not  converge  to  the  true  parameter  values.  We  now  employ  the  diagonal 
pertubation  technique  in  a  manner  similar  to  that  of  Chapter  III  which  improves  the 
condition  of  R  and  permits  a  modified  version  of  (4.8)  to  converge  as  the  number  of 
data  points  approaches  the  minimum  required. 

If  we  let  (7y  and  Cx  be  scalar  multipliers  and  the  Toeplitz  approximation 
splitting  of  Ryy  and  Rxx  be  defined  as  Ryy  —  Tyy^  -f  Syy  and  R^x  =  Txxo  + 
respectively,  then  we  can  derive  a  modified  version  of  (4.8)  as  follows: 

R0  =T  (4.10) 


adding  a  diagonal  matrix  to  Ryy  and  Rxx  ijives, 


Ryy  +  (Tyl 

:  Rj/x 

b' 

ry' 

+ 

Vyb' 

Rxy 

•  Rxx 

a 

<Txa_ 

and  substituting  for  Ryy  and  Rxx  we  have 


(4.11) 


TyVo  +  Syy  +  CTy/ 

Ryx 

h' 

ry' 

+ 

cTyb' 

Rxy 

Txxo  +  Sxx  +  CTxI  _ 

a 

<7xa 

(4,12) 


Now,  letting  Tyy  =  Tyy^  +  Cyl  and  Txx  =  Txxo  +  and  multiplying  gives  the 
following  coupled  equations: 


{Tyy  +  5vy)b  +  Ryx^  =  +  CTyb 

{Txx  +  5ir)a  +  Rxyh  =  Tx  +  cTxa  (4-13) 


and  rearranging  these  gives 
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(4.14) 


Tyyh  =  Tj,  +  (Tyh  -  Ry^a  -  Syyh 

I'xx^  -  4”  .^xy^  Sxx^‘ 

Substituting  for  Syy  =  Ryy  +  (Tyl  —  Tyy  and  Sn  =  Rxx  +  cTxI  —  T^x  we  have 

*  ry"|"^yb  Ry^a  {^Ryy  “j”  ^yf  ^yy)^ 

Tx^a  =  Tx  +  C7xa  -  i^xyb  -  (i?xx  +  Cx-f  -  ?’xx)a,  (4-15) 

and  now  we  see  the  terms  £rj,b  and  a^a  cancel  out: 

Tyyb  —  ry"[“^yb  ^yb  Ryx^  ("^^yy  ^yy)^ 

Txxa  =  Tx  +  <Txa  -  (Txa  -  i?xyb  -  (i?xx  -  2^xx)a.  (4.16) 

After  simplifying  to  isolate  the  parameter  vectors,  we  have 

b  =  T-^^y  -  T~^{Ryxa  +  {Ryy  "  rj,j,)b) 

a  =  -  T^x^Rxyh  +  {Rxx  -  Txx)^)-  (4.17) 

Again,  letting  T^^Vy  =  bo  and  T~^rx  =  ao  gives 

b  =  bo  +  b-r~^(i?j,xa  +  i2yjb) 

a  =  ao  +  a  -  7’;;^(iZxyb  +  i^ixa)  (4.18) 

amd  the  finad  form  of  the  Toeplitz  approximation  with  diagonad  pertubation  iterative 
algorithm  is: 

=  bo  +  b(*^  -  r-'(i2„xa^^)  +  i?y,b('')) 

=  ao  +  aW-r;;'(i?xyb(*+')  +  i?xxa(*')).  (4.19) 

The  algorithm  of  (4.19)  hais  exactly  the  saune  form  as  (4.8).  Unfortunately,  adding 
a  diagonal  value  dramatically  slows  down  the  rate  of  convergence;  however,  the  new 
adgorithm  (4.19)  was  observed  to  converge  where  the  algorithm  (4.8)  fails.  This 
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Figure  4.2:  Parameter  tracks  and  general  rate  of  convergence  for  the  IIR 
Toeplitz  approximation  with  diagonal  pertubation  algorithm. 

algorithm  takes  a  large  number  (i.e.,  several  hundreds  or  thousands)  of  iterations  just 
to  reach  the  vicinity  of  the  true  parameter  values,  and  still  requires  a  large  number  of 
further  iterations  to  converge  to  the  correct  solution.  We  did  observe  that  appropriate 
values  for  <Ty  and  <7*  are  inversely  proportional  to  the  condition  number  of  Ryy  and 
Rxxt  respectively.  Unfortunately,  even  the  optimal  performance  of  this  algorithm  is 
unsatisfactory  (see  Figure  4.2). 

3.  Partitioned  Toeplitz  Approximation  Algorithm 

In  this  section  we  block  partition  Ryy  and  Rxx  in  an  attempt  to  improve  the 
performance  of  (4.8).  As  before,  the  number  of  row  partitions  must  equ2il  the  number 
of  colunm  partitions.  Similar  to  previous  developments  we  now  derive  a  partitioned 
version  of  (4.8).  Beginning  with 

Re  =  r,  (4.20) 


0  5000  10000 

Iterations,  k 
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partitioning  twice  gives: 


(■^^11  I  ^yi2  \  f  -^*11  I  -^*12  '\  ^vi 

^y21  1  •^»22  /  V  ^121  1  ■^X22  /  ^2  ^^2 

/  -Raryji  |  ^xyi2  \  /  ^xxli  I  -^1x12  \  ^^xl 

\  ■^XV21  I  •^11/22  /  \  ■^^21  I  Rxx22  /  .  .  ^2  _  _  ^^2 

Splitting  the  diagonaJ  blocks  into  a  Toeplitz  and  a  residual  matrix  we  have 

f  (  '^yy\\  d"  ‘^yvii  I  ^yi2  \  (  -^^ii  I  -^*12  \ 


J^xyii  I  ^xyi2 


Rxy2l  I  •^XV22 


■yvii  ~  “^yvii 
^y2i 


'^yy22  *^^22 


•^IV21  I  ^^y22 


(  AfXll 

1  -^112 

\  Ryx2\ 

1  -^X22 

Txx\\  d"  *Sixll 

1  R: 

Rxx21 

1  Txx22  ' 

(4.21) 


Multiplying  gives  four  coupled  equations: 


(4.22) 


(^VVll  ‘^yVll)^l  "h  -^^12^2  d"  -^xii^l  d"  .f^xi2^2  —  Tyi 

(^yy22  ‘^yy22)^2  •^^yv2i^^  •l^x2i^i  d"  Ryx22^2  —  ^112 

(^xxll  d- ■S'r*ii)ai  +  iZiy^lbi  + /?xvi2^2  d- i2xxl2®2  =  fxl 

(2*x22  d"  ‘^xx22)^2  d"  •^2x^21  bi  d"  ■^xy22^^  ■^^x21®l  ~  ^x2  (4.23) 

and  rearranging  yields 
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~  Ryxii^l  Ryxi2^2 

^VV22^^  “  ^1/2  •^^y21^1  •^^*22®'^  ‘^yy22^^ 

^rxll^l  —  -Rrviibi  Rxyy^2  •^®xl2^2  •S'xxH^l 

1'xx22^3  ~  ^x2  •^ry21^1  -^13/22^2  ^xi21®l  ‘^ix22®’2'  (‘^•24) 

Now,  isolating  the  parameter  vectors  and  substituting  for  Syy^^  =  Hyy^^  —  Tyy^^,  etc., 

bi  =  Tyy  —  Tyy  j j ( -^ V 1 2 ^ ^  ”  ^^yx j j "  -^112^^2  ~  (-^yii  ~  ^yyii)^i) 

^>2  =  T'yy  2'^y2  ~  '^W  22('^S'21^^  ~  ■^121^1  ~  -^122^2  “  {Ryy22  ~  ^yy22)^2) 

Hi  =  Tj-x  11^x1  “  ^xx  llC-^^ryil^l  “  •Rxyi2^2  ~  Rxx\2^2  ~  (-^rxll  “  Jxill)3-l) 

^2  “  Txx  22*^®2  ~  ^ii  22('^3^y21^’-  ~  ■^®y22^2  ~  .Rix21®l  ""  (■^xi22  ~  ^xx22)^2)  (4.25) 

and  letting  bio  =  ^yy^n'^vi’  ^  initial  estimate  of  the  parameter  vectors,  we 
have 

bi  bio  "h  ^yy  ii(-^^yyi2^2  ~  -Ryxn^l  ~  Ryxi2^^  -^^^yyii^i^) 

b2  =  b2o  4"  b2  —  Tyy  22(-^^^21^1  ~  ■^X2i^l  -^X22®2  •^'^22^2) 

ni  —  Sio  4"  Si  T^x  ii(-^iyil^l  •^xyi2^2  -^^112^2  -^^xll^l) 

S2  =  S2o  4-  a2  —  Txx^ 22(-^^21^1  ~  ■^y22^2  ~  ■^x2lSi  —  Rxx22^2)'  (4-26) 

The  final  form  of  the  partitioned  fixed  data  Toeplitz  approximation  iterative  algorithm 
is: 

bS*+‘>  =  b.„  +  b[‘>  -  r-‘„(A„„b5‘>  - 

b‘‘+'>  =  b,„  +  b'"'  -  r-‘„(il„,.b<‘«'  -  -  /J„„b<‘>) 

al*+'>  =  aio  +  af  -  -  i?„„b5*+'>  -  iJ.„,af  -  fi.,„al*') 

a5**”  =  ajo  +  ai‘>  -  T,->„(A,„bS‘"'  -  fl„„b<‘+‘l  -  ie,x„a<‘-"'>  -  iJ.„,aS‘>). 

(4.27) 

This  partitioned  Toeplitz  approximation  algorithm  w«is  observed  to  con¬ 
verge  to  the  true  solution  where  the  algorithm  (4.8)  fails,  however  it  hais  the  same 
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Figure  4.3:  Parameter  tracks  and  general  rate  of  convergence  for  the  HR 
partitioned  Toeplitz  approximation  algorithm. 

problem  as  the  diagonal  pertubation  technique.  It  takes  an  excessive  number  of  itera¬ 
tions  for  the  algorithm  to  approach  the  vicinity  of  the  true  parameter  values  and  even 
longer  to  converge  exactly  (see  Figure  4.3).  Since  this  is  not  acceptable  performance 
(even  marginally),  we  therefore  turn  to  data-adaptive  techniques  in  order  to  achieve 
the  desired  result. 

B.  DATA-ADAPTIVE  ALGORITHMS 
1.  Toeplitz  Approximation  Algorithm 

Extending  the  development  from  equations  (4.1)-(4.6)  to  the  time- varying 
case  gives  the  desired  data-adaptive  IIR  Toeplitz  approximation  algorithm: 

bl‘«'  =  b(‘>  + 

a(*+»  =  a(‘)+r„;‘(r„-ft,„b<‘)-iJ„„aW)  (4.28) 

where  the  subscript  n  denotes  time  variation.  The  neat  and  clean  appearance  of 
this  algorithm  is  deceptive.  Each  of  the  terms  with  a  subscript  n  must  be  updated 
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during  each  iteration.  Consequently,  this  algorithm  is  somewhat  computationally 
intensive.  The  entire  algorithm  is  coded  in  matlab  and  is  listed  in  Appendix  D.  It 
uses  the  matrix  inversion  lemma  to  compute  the  Tyy^  and  Txxn  updates.  Updating 
the  rest  of  the  time  varying  terms  is  similar  to  the  process  shown  in  (3.37)  and 
(3.38).  We  observed  that  the  algorithm  gave  better  results  for  over-modeling  than 
exact  modeling.  An  example  of  over-modeling  is  where  the  true  system  output  was 
generated  from  a  4*^-order  system  {M  =  N  =  4),  but  we  are  using  a  6‘^-order  model 
(M  =  iV  =  6)  in  the  algorithm.  The  results  for  this  case  are  shown  in  Figure  4.4. 
Here  we  see  that  the  parameter  values  converge  in  about  15  to  20  iterations.  To 
see  whether  this  6‘^-order  model  faithfully  represents  the  true  4*^-order  system,  we 
must  look  at  the  frequency  response  of  the  model  versus  the  true  system.  Figure 
4.5  shows  that  by  about  20  iterations  the  6*^-order  model  matches  the  frequency 
response  of  the  true  system.  This  performance  is  quite  reasonable,  and  Appendix  C 
contains  parameter  track  and  frequency  response  plots  showing  its  performance  under 
a  variety  of  conditions.  Another  important  point  about  this  algorithm  is  that  it  does 
not  use  the  gradient  terms  that  appe2ir  in  (2.21).  In  fact,  when  the  gradient  terms 
were  incorporated  in  the  algorithm,  the  performance  was  observed  to  be  worse. 

2.  Recursive  Least  Squares  (RLS)  Algorithm 

Unfortunately,  there  is  not  a  well-defined  algorithm  that  will  serve  as  a 
benchmark  against  which  to  test  the  data-adaptive  Toeplitz  approximation  algorithm. 
In  this  section  we  use,  without  derivation,  the  IIR  RLS  algorithm  defined  by[Ref.  5]: 

0(fc+i)  =  gW  +  aRzUMn)  (4.29) 

where  6  =  [b^  a^]^  is  the  parameter  vector,  is  the  gradient  term 

vector,  and  a  is  the  forgetting  factor.  There  are  some  real  difficulties  in  keeping  this 
algorithm  stable  [Ref.  5],  and  it  performs  best  when  used  to  exactly  model  a  system’s 
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Figure  4.4:  Parameter  tracks  and  general  rate  of  convergence  for  the  IIR 
data>adaptive  Toeplitz  approximation  algorithm. 
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Figure  4.5:  Frequency  response  of  the  Toeplitz  approximation  algorithm 
compared  to  the  true  system  frequency  response. 
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Figure  4.6:  Parameter  tracks  and  general  rate  of  convergence  for  the  IIR 
RLS  algorithm. 

parameters.  The  performance  of  this  algorithm,  when  it  is  stable,  is  shown  in  Figure 
4.6.  However,  for  general  system  modeling  where  the  exact  order  of  the  system  is 
unknown,  IIR  RLS  is  unsatisfactory.  See  Figures  C.IO  thru  C.14  of  Appendix  C. 

Comparing  the  performance  of  IIR  RLS  to  that  of  data-adaptive  IIR  Toeplitz 
approximation,  we  observe  that  the  increaised  computations  required  in  the  Toeplitz 
approximation  edgorithm  give  improved  performance.  The  stability  of  the  Toeplitz 
approximation  algorithm  is  cleaxly  evident  when  we  consider  the  poles  and  zeros  of 
the  model  system  as  it  converges.  In  Figure  4.7  the  poles  (denoted  by  x)  amd  zeros 
(denoted  by  o)  move  inside  the  unit  circle  ais  the  number  of  iterations  goes  from  9  to 
20.  The  fact  that  the  poles  and  zeros  stop  moving  as  the  algorithm  continues  to  iter¬ 
ate  cleaurly  demonstrates  that  it  is  stable.  Also  observe  that  the  two  additional  poles 
and  zeros  overlap  at  iteration  k  =  50,  which  effectively  cancels  their  contribution  to 
the  frequency  response. 

Similair  performance  is  observed  in  the  presence  of  noise,  and  the  Toeplitz 
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Figure  4.8:  Parameter  track  and  frequency  response  for  the  data-adaptive 
HR  Toeplitz  approximation  algorithm  with  M  =  N  =  20  and  10  dB  additive 
noise. 

approximation  algorithm  converges  faster  to  the  true  frequency  response  when  higher 
order  model  systems  are  computed.  Notice  how  the  true  frequency  response  emerges 
as  the  algorithm  begins  to  converge  for  a  M  =  iV  =  20  model  system  in  the  presence 
of  10  dB  additive  noise  (see  Figure  4.9). 
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V.  CONCLUSIONS 

In  this  thesis  we  applied  the  four  concepts  identified  in  the  Introduction  to 
both  FIR  and  HR  system  parameter  estimation.  Of  these,  diagonal  pertubation  and 
multiple  paxtitioning  are  used  to  extend  the  basic  fixed  data  algorithms  of  [Ref.  1] 
for  FIR  systems  and  of  [Ref.  2]  for  HR  systems.  We  then  used  all  four  concepts  to 
create  new  data-adaptive  algorithms  which  extend  the  fixed  data  algorithms  into  the 
data-adaptive  domain.  In  the  following  sections  we  draw  general  conclusions  about 
the  performance  of  these  algorithms. 

A.  FIR  ALGORITHMS 

Improving  the  condition  of  the  covariajice  matrix  with  increcised  diagonal  dom¬ 
inance  produced  improved  performance  over  the  original  fixed  data  FIR  Toeplitz 
approximation  algorithm.  The  algorithm  converges  fairly  rapidly  while  using  the 
minimum  eunount  of  data.  Figures  A.l  thru  A.5  of  Appendix  A  show  the  general  per¬ 
formance  of  this  aJgorithm  in  noise,  ais  well  as  for  over-modeling  and  under-modeling. 
The  major  drawback  of  this  algorithm  is  the  problem  of  determining  a  suitable  vailue 
of  the  diagonal  pertubation  parameter  a  that  will  give  optimad  performance. 

Partitioning  the  covariance  matrix  also  produced  improved  performance  of  the 
originaJ  fixed  data  FIR  Toeplitz  approximation  algorithm.  Using  the  minimum  amoimt 
of  data,  the  partitioned  Toeplitz  approximation  algorithm  converges  almost  as  rapidly 
as  the  optimal  case  of  the  diagonal  pertubation  Toeplitz  approximation  zilgorithm. 
Figures  A. 6  thru  A. 9  of  Appendix  A  show  the  generad  performance  of  this  algorithm 
in  noise,  <is  well  as  for  over-modeling  and  under-modeling.  Significant  features  of  this 
algorithm  are  that  it  does  not  require  choosing  am  optimizing  pwameter,  and  the 
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paxtitioning  reduces  the  computational  intensity  of  the  algorithm  since  the  matricies 
involved  zu'e  smaller. 

Extending  the  original  fixed  data  FIR  Toeplitz  approximation  algorithm  to  the 
data-adaptive  case  gives  an  algorithm  that  performs  quite  well,  equal  to  RLS,  however 
it  is  more  computationally  intensive  than  RLS.  Establishing  this  data-adaptive  FIR 
algorithm  is  important  because  it  demonstrates  the  potential  for  a  data-adaptive  HR 
version  to  work  as  well.  Figures  A.  10  thru  A.  17  of  Appendix  A  show  the  general 
performance  of  this  algorithm  and  RLS  under  different  conditions  of  noise,  cis  well  as 
for  over-modeling  and  under-modeling. 

B.  IIR  ALGORITHMS 

Applying  the  same  diagonal  pertubation  and  partitioning  techniques  used  on 
the  fixed  data  FIR  Toeplitz  approximation  algorithm  to  the  IIR  algorithm  did  not 
produce  the  same  improvement.  Although  these  techniques  did  allow  the  algorithm 
to  converge  with  the  desired  minimum  amount  of  data,  the  rate  of  convergence  is 
unacceptable.  Figures  B.l  thru  B.9  of  Appendix  B  show  the  general  performance 
of  these  algorithms  under  different  conditions  of  data  md  noise,  eis  well  as  for  over¬ 
modeling  and  under-modeling. 

The  data-adaptive  IIR  Toeplitz  approximation  algorithm  succeeds  in  parameter 
estimation  where  the  corresponding  IIR  RLS  algorithms  feiil.  It  performs  well  under  a 
wide  range  of  noise  and  over-modeling  conditions.  We  observed  that  even  in  moderate 
noise  (i.e.,  10  to  20  dB)  over-modeling  with  this  adgorithm  can  be  used  to  determine 
the  true  frequency  response  of  the  system  being  modeled.  Figures  C.l  thru  C.14  of 
Appendix  C  show  how  the  2dgorithms  perform  under  different  conditions  of  noise, 
as  well  as  over-modeling  and  under-modeling.  A  significant  consideration  with  this 
algorithm  is  the  large  amount  of  computations  required  for  each  iteration. 
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C.  FUTURE  WORK 


Here  we  have  reported  general  observations  rather  than  a  rigorous  analytical 
development  on  the  performance  of  the  algorithms  presented  in  the  thesis.  As  future 
work,  the  results  and  conclusions  from  this  work  should  be  put  on  solid  analytical 
footing.  Also,  the  algorithms  presented  in  this  thesis  are  developed  for  one-dimensonal 
real  data  only.  The  algorithms  may  be  extended  to  incorporate  complex  and  multi¬ 
dimensional  data  and  related  applications. 
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APPENDIX  A:  FIR  SYSTEM  SIMULATIONS 


This  appendix  contains  plots  corresponding  to  the  FIR  algorithms  developed  in 
Chapter  III.  Figure  A.l  shows  the  performance  of  the  unmodified  Toeplitz  approxi¬ 
mation  algorithm  (3.15).  This  plot  also  serves  as  a  reference  of  the  true  system  (3.8) 
frequency  response.  The  remainder  of  the  plc^s,  Figures  A.2-A.17,  are  in  the  follow¬ 
ing  general  order:  fixed  data  Toeplitz  approximation  with  diagonal  pertubation,  fixed 
data  partitioned  Toeplitz  approximation,  data-adaptive  Toeplitz  approximation,  and 
finally  FIR  RLS. 
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Figure  A.l:  The  basic  fixed  data  FIR  Toeplitz  approximation  algorithm. 
This  plot  shows  the  parameter  tracks  and  frequency  response  of  the  baisic  Toeplitz 
approximation  algorithm  (3.15).  The  correlation  matrix  was  formed  with  P  =  20 
data  points,  and  there  wets  no  additive  noise  in  the  output.  The  algorithm  converges 
rapidly  under  these  conditions,  as  can  be  seen  by  the  flat  p2ir2uiieter  tracks.  Also 
notice  that  the  frequency  response  changes  little  as  the  number  of  iterations  increase. 
The  frequency  response  plot  for  k  =  20  iterations  represents  the  frequency  response 
of  the  true  system. 
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Figure  A. 2:  Parameter  tracks  and  frequency  response  of  the  fixed  data 
Toeplitz  approximation  with  diagonal  pertubation  algorithm. 

Observe  that  even  though  the  minimum  number  of  data  points  (P  =  11)  were  used  to 
form  the  correlation  matrix,  the  algorithm  converges  rather  quickly.  The  frequency 
response  plot  for  fc  =  30  iterations  appeaurs  identical  to  the  response  of  the  true  system 
shown  in  Figure  A.l.  No  noise  was  added  to  the  output,  the  pertubation  parameter 
was  <T  =  1,  2md  the  true  system  was  exactly  modeled  with  a  10‘^-order  model. 
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Figure  A. 3:  Under-modeled  example  of  the  fixed  data  Toeplitz  approxi¬ 
mation  with  diagonal  pertubation  algorithm. 

The  10‘^-order  true  system  was  modeled  with  £in  8*^-order  model.  The  algorithm  does 
a  fairly  good  job  preserving  the  frequency  components  of  the  true  system,  however  it 
takes  a  much  greater  number  of  iterations  to  achieve  this  result.  No  noise  was  added 
to  the  output,  and  the  pertubation  parameter  was  (7=1. 
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Figure  A. 4:  Over-modeled  example  of  the  fixed  data  Toeplitz  approxima¬ 
tion  with  diagonal  pertubation  algorithm. 

The  10*^-order  true  system  was  modeled  with  a  14‘*-order  model.  The  eilgorithm 
converges  rapidly,  said  by  A:  =  30  iterations  it  matches  the  true  system  frequency 
response.  It  is  important  to  note  that  this  result  was  achieved  with  the  minimum 
[P  =  15)  amount  of  data  required.  No  noise  was  added  to  the  output,  <ind  the 
pertubation  parsuneter  was  <t  =  3. 
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Figure  A. 5;  Over-modeled  example  with  10  dB  of  additive  noise  for  the 
fixed  data  Toeplitz  approximation  with  diagonal  pertubation  algorithm. 
The  lO^^-order  true  system  was  modeled  with  a  20‘'‘-order  model.  The  additive  white 
noise  tends  to  slow  down  the  algorithm,  and  distorts  the  frequency  response  it  is  able 
to  achieve.  The  pertubation  parameter  was  (7  =  3,  and  10  dB  of  white  noise  was 
added  to  the  output. 
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Figure  A.6:  Parameter  tracks  and  frequency  response  of  the  fixed  data 
partitioned  Toeplitz  approximation  algorithm. 

Observe  that  even  though  the  minimum  number  of  data  points  (P  =  11)  were  used  to 
form  the  correlation  matrix,  the  £ilgorithm  converges  rather  quickly.  The  frequency 
response  plot  for  k  =  30  iterations  appears  identical  to  the  response  of  the  true  system 
shown  in  Figure  A.l.  No  noise  was  added  to  the  output,  and  the  true  system  was 
exactly  modeled  with  a  10‘'‘-order  model. 
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Figure  A. 7:  Under-modeled  example  of  the  fixed  data  partitioned  Toeplitz 
approximation  algorithm. 

The  10*^-order  true  system  was  modeled  with  an  8*^-order  model.  The  algorithm  does 
a  faurly  good  job  preserving  the  frequency  components  of  the  true  system,  however  it 
takes  a  much  greater  number  of  iterations  to  achieve  this  result.  No  noise  was  added 
to  the  output. 
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Figure  A. 8:  Over-modeled  example  of  the  fixed  data  partitioned  Toeplitz 
approximation  algorithm. 

The  10‘^-order  true  system  was  modeled  with  a  14*^-order  model.  The  «ilgorithm 
converges  rapidly,  and  by  t  =  30  iterations  it  matches  the  true  system  frequency 
response.  It  is  important  to  note  that  this  result  was  achieved  with  the  minimum 
{P  =  15)  amount  of  data  required.  No  noise  was  added  to  the  output. 
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Figure  A. 9:  Over-modeled  example  with  10  dB  of  additive  noise  for  the 
fixed  data  partitioned  Toeplitz  approximation  algorithm. 

The  10‘^-order  true  system  was  modeled  with  a  20*^-order  model.  The  additive  white 
noise  tends  to  slow  down  the  algorithm,  and  distorts  the  frequency  response  it  is  able 
to  achieve.  There  was  10  dB  of  white  noise  added  to  the  output. 
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Figure  A.  10:  Parameter  tracks  and  frequency  response  of  the 

data-adaptive  Toeplitz  approximation  algorithm. 

Observe  that  the  algorithm  converges  to  the  true  solution  and  frequency  response  in 
ten  iterations.  In  this  example  no  noise  was  added  to  the  output,  and  the  performance 
above  is  for  an  exact  modeling  of  the  10‘'*-order  true  system.  Comparing  this  with 
the  RLS  results  shown  in  Figure  A.  14  reveals  that  they  behave  exactly  the  same. 
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Figure  A.  11;  Under-modeled  example  of  the  data-adaptive  Toeplitz  ap¬ 
proximation  algorithm. 

The  10**-order  true  system  was  modeled  with  an  8‘*-order  model.  The  algorithm 
tries  to  preserve  the  frequency  components  of  the  true  system,  however  it  taJses  a 
much  greater  number  of  iterations  to  achieve  this  result.  No  noise  was  added  to  the 
output. 
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Figure  A. 12:  Over-modeled  example  of  the  data-adaptive  Toeplitz  approx¬ 
imation  algorithm. 

The  10‘*-order  true  system  was  modeled  with  a  14‘^-order  model.  The  algorithm 
converges  rapidly,  and  for  iterations  «  =  10  and  higer  it  matches  the  true  system 
frequency  response.  It  is  important  to  note  that  aJl  of  the  additional  pareimeters 
beyond  the  ten  required  were  driven  to  zero  by  the  algorithm.  No  noise  was  added 
to  the  output. 
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Figure  A.  13;  Over-modeled  example  with  10  dB  of  additive  noise  for  the 
data-adaptive  Toeplitz  approximation  algorithm. 

The  10‘^-order  true  system  was  modeled  with  a  20*^-order  model  The  additive 
white  noise  tends  to  slow  down  the  algorithm,  but  only  slightly  distorts  the  frequency 
response  it  is  able  to  aw:hieve.  There  was  10  dB  of  white  noise  added  to  the  output. 
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Figure  A. 14:  Parameter  tracks  and  frequency  response  of  the  FIR  RLS 
algorithm. 

Observe  that  the  algorithm  converges  to  the  true  solution  and  frequency  response  in 
ten  iterations.  In  this  example  no  noise  was  added  to  the  output,  and  the  performance 
above  is  for  tin  exact  modeling  of  the  10‘^-order  true  system.  Comparing  this  with  the 
data-adaptive  Toeplitz  approximation  algorithm  results  shown  in  Figure  A. 10  reveals 
that  they  behave  exactly  the  same. 
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Figure  A. 15:  Under-modeled  example  of  the  FIR  RLS  algorithm. 

The  10‘^-order  true  system  was  modeled  with  an  8*^-order  model.  The  algorithm 
tries  to  preserve  the  frequency  components  of  the  true  system,  however  it  takes  a 
much  greater  number  of  iterations  to  achieve  this  result.  No  noise  wcis  added  to  the 
output. 
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Figure  A.  16:  Over-modeled  example  of  the  FIR  RLS  algorithm. 

The  10‘^-order  true  system  was  modeled  with  a  14‘^-order  model.  The  algorithm 
converges  rapidly,  and  by  iteration  A:  =  10  and  higher  it  matches  the  true  system 
frequency  response.  It  is  important  to  note  that  all  of  the  additional  parameters 
beyond  the  ten  required  were  driven  to  zero  by  the  algorithm.  No  noise  was  added 
to  the  output. 
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Figure  A. 17:  Over-modeled  example  with  10  dB  of  additive  noise  for  the 
data-adaptive  Toeplitz  approximation  algorithm. 

The  10*^-order  true  system  wm  modeled  with  a  20‘*-order  model.  The  additive 
white  noise  tends  to  slow  down  the  algorithm,  but  only  slightly  distorts  the  frequency 
response  it  is  able  to  achieve.  There  was  10  dB  of  white  noise  added  to  the  output. 
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APPENDIX  B:  FIXED  DATA  HR  SYSTEM 

SIMULATIONS 

This  appendix  contains  plots  corresponding  to  the  fixed  data  HR  algorithms  de¬ 
veloped  in  Chapter  IV.  Figure  B.I  shows  the  performance  of  the  unmodified  Toeplitz 
approximation  algorithm  (4.8).  This  plot  also  serves  as  a  reference  of  the  true  sys¬ 
tem  (4.9)  frequency  response.  The  remainder  of  the  plots,  Figures  B.2-B.9  illustrate 
the  convergence  performance  of  the  fixed  data  Toeplitz  approximation  with  diagonal 
pertubation  algorithm  followed  by  the  fixed  data  partitioned  Toeplitz  approximation 
algorithm. 
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Figure  B.l:  The  basic  fixed  data  IIR  Toeplitz  approximation  algorithm. 
This  plot  shows  the  parameter  tracks  and  frequency  response  of  the  basic  Toeplitz 
approximation  algorithm  (4.8).  The  correlation  matrix  wets  formed  with  P  =  100 
data  points,  and  there  W2is  no  additive  noise  in  the  output.  The  algorithm  converges 
rapidly  under  these  conditions,  as  can  be  seen  by  the  flat  parameter  tracks.  Also 
notice  that  the  frequency  response  changes  little  «is  the  number  of  iterations  increase. 
The  frequency  response  plot  for  k  —  ZO  iterations  represents  the  frequency  response 
of  the  true  system. 
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Figure  B.2:  Parameter  tracks  and  frequency  response  of  the  fixed  data 
Toeplitz  approximation  with  diagonal  pertubation  algorithm. 

Observe  that  when  the  minimum  number  of  data  points  (P  =  9)  were  used  to  form 
the  correlation  matrix,  the  algorithm  converges  very  slowly.  The  frequency  response 
plot  for  k  =  50  iterations  is  not  very  close  to  the  response  of  the  true  system  shown 
in  Figure  B.l.  No  noise  was  added  to  the  output,  the  pertubation  parameters  were 
(T,  =  1  and  (Ty  =  10,  and  the  true  system  was  exactly  modeled  with  a  4‘^-order  model. 
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Figure  B.3:  Under-modeled  example  of  the  fixed  data  Toeplitz  approxi¬ 
mation  with  diagonal  pertubation  algorithm. 

The  4'*-order  true  system  wm  modeled  with  an  3’''^-order  model.  The  algorithm  does 
a  poor  job  preserving  the  frequency  components  of  the  true  system,  and  it  takes  a 
much  greater  number  of  iterations  to  achieve  this  result.  No  noise  was  added  to  the 
output,  and  the  pertubation  parameters  were  (Ti  =  4  ^lnd  <Ty  =  50. 
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Figure  B.4:  Over-modeled  example  of  the  fixed  data  Toeplitz  approxima¬ 
tion  with  diagonal  pertubation  algorithm. 

The  4‘^-order  true  system  was  modeled  with  a  10‘^-order  model.  The  algorithm 
converges  more  rapidly  than  the  previous  two  examples,  and  by  fc  =  100  iterations  it 
nearly  matches  the  true  system  frequency  response.  It  is  important  to  note  that  this 
result  was  achieved  with  the  minimum  (P  =  21)  amount  of  data  required.  No  noise 
was  added  to  the  output,  and  the  pertubation  parameters  were  —  Z  and  o-y  =  100. 
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Figure  B.5;  Over-modeled  example  with  10  dB  of  additive  noise  for  the 
fixed  data  Toeplitz  approximation  with  diagonal  pertubation  algorithm. 
The  4‘*-order  true  system  was  modeled  with  a  20*^-order  model.  The  additive  white 
noise  tends  to  slow  down  the  algorithm,  but  only  slightly  distorts  the  frequency 
response  it  is  able  to  achieve.  The  pertubation  pareimeters  were  cTj.  =  5  and  ay  =  100, 
and  10  dB  of  white  noise  was  added  to  the  output. 
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Figure  B.6;  Parameter  tracks  and  frequency  response  of  the  fixed  data 
partitioned  Toeplitz  approximation  algorithm. 

Observe  that  when  the  minimum  number  of  data  points  (P  =  9)  were  used  to  form 
the  correlation  matrix,  the  algorithm  converges  very  quickly.  The  frequency  response 
plot  for  fc  =  50  iterations  is  not  very  close  to  the  response  of  the  true  system  shown 
in  Figure  B.l.  No  noise  was  added  to  the  output,  aind  the  true  system  was  exactly 
modeled  with  a  4‘'‘-order  model. 
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Figure  B.7:  Under-modeled  example  of  the  fixed  data  partitioned  Toeplitz 
approximation  algorithm. 

The  4‘^-order  true  system  was  modeled  with  cin  3'’‘^-order  model.  The  algorithm  does 
a  poor  job  preserving  the  frequency  components  of  the  true  system,  and  it  takes  a 
much  greater  number  of  iterations  to  achieve  this  result.  No  noise  wcis  added  to  the 
output. 
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Figure  B.8:  Over-modeled  example  of  the  fixed  data  partitioned  Toeplitz 
approximation  algorithm. 

The  4*^-order  true  system  was  modeled  with  a  10‘*-order  model.  The  algorithm 
converges  rapidly,  smd  by  t  =  20  iterations  it  matches  the  true  system  frequency 
response.  It  is  important  to  note  that  this  result  wcis  achieved  with  the  minimum 
(P  =  21)  eunount  of  data  required.  No  noise  was  added  to  the  output. 
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Figure  B.9:  Over-modeled  example  with  10  dB  of  additive  noise  for  the 
fixed  data  partitioned  Toeplitz  approximation  algorithm. 

The  4*^-order  true  system  was  modeled  with  a  20*^-order  model.  The  additive  white 
noise  tends  to  slow  down  the  algorithm,  but  only  slightly  distorts  the  frequency 
response  it  is  able  to  achieve.  There  was  10  dB  of  white  noise  added  to  the  output. 
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APPENDIX  C:  DATA-ADAPTIVE  HR 
SYSTEM  SIMULATIONS 


This  appendix  contains  plots  corresponding  to  the  data-adaptive  HR  algorithms 
developed  in  Chapter  IV.  Figure  C.l  shows  the  performance  of  the  unmodified  fixed 
data  Toeplitz  approximation  algorithm  (4.8).  This  plot  also  serves  as  a  reference  of 
the  true  system  (4.9)  frequency  response.  The  remainder  of  the  plots,  Figures  C.2- 
C.4,  illustrate  the  performance  of  the  data-adaptive  Toeplitz  approximation  algorithm 
followed  by  the  IIR  RLS  algorithm. 
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Figure  C.l:  The  basic  fixed  data  HR  Toeplitz  approximation  algorithm. 
This  plot  shows  the  parameter  tracks  and  frequency  response  of  the  beisic  Toeplitz 
approximation  algorithm  (4.8).  The  correlation  matrix  was  formed  with  P  =  100 
data  points,  and  there  wjis  no  additive  noise  in  the  output.  The  cilgorithm  converges 
rapidly  under  these  conditions,  as  can  be  seen  by  the  flat  parameter  tracks.  Also 
notice  mat  the  frequency  response  changes  little  eis  the  number  of  iterations  increase. 
The  frequency  response  plot  for  k  =  30  iterations  represents  the  frequency  response 
of  the  true  system. 
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Figure  C.2;  Parameter  tracks  and  frequency  response  of  the  data-adaptive 
Toeplitz  approximation  algorithm. 

Observe  that  the  algorithm  converges  to  the  true  solution  and  frequency  response 
in  about  60  iterations.  In  this  example  no  noise  was  added  to  the  output,  and  the 
performance  above  is  for  an  exact  modeling  of  the  4*'‘-order  true  system.  Comparing 
this  with  the  RLS  results  shown  in  Figure  B.18  reveals  that  they  do  not  behave  the 
same  as  in  the  FIR  case. 
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Figure  C.3:  Pole-zero  plots  for  the  data-adaptive  Toeplitz  approximation 
algorithm. 

Observe  that  as  the  algorithm  converges  to  the  true  solution,  the  poles  (x’s)  move 
inside  the  unit  circle.  After  about  50  iterations  they  stop  moving,  and  their  location 
corresponds  to  the  frequency  response  shown  in  Figure  B.IO.  Finally,  we  see  that 
there  is  almost  no  movement  as  the  algorithm  iterates  from  ifc  =  70  to  100  iterations. 
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Figure  C.4:  Under- modeled  example  of  the  data-adaptive  Toeplitz  approx¬ 
imation  algorithm. 

The  4**-order  true  system  was  modeled  with  a  3’’'^-order  model.  The  algorithm  tries 
to  preserve  the  frequency  components  of  the  true  system,  but  with  only  three  poles 
aind  zeros  it  is  unable  to  preserve  both  of  the  m{un  components.  No  noise  was  added 
to  the  output. 
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Figure  C.5:  Pole-zero  plots  for  the  under-modeled  data-adaptive  Toeplitz 
approximation  algorithm. 

Observe  that  as  the  algorithm  converges,  the  poles  (x’s)  still  move  inside  the  unit 
circle.  After  about  50  iterations  their  location  corresponds  to  the  frequency  response 
shown  in  Figure  B.12.  Finally,  we  see  that  the  pole-zero  pair  on  the  real  axis  tend  to 
cancel  each  other,  which  explains  why  there  is  only  one  primary  frequency  component 
in  Figure  B.12. 
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Figure  C.6:  Over-modeled  example  of  the  data-adaptive  Toeplitz  approx¬ 
imation  algorithm. 

The  4**-order  true  system  was  modeled  with  a  10‘^-order  model.  The  algorithm 
converges  rapidly,  and  by  it  =  50  iterations  it  matches  the  true  system  frequency 
response.  It  is  important  to  note  that  all  of  the  additional  pole  zero  pairs  beyond  the 
four  required  were  driven  to  cancel  each  other  by  the  algorithm.  No  noise  was  added 
to  the  output. 
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Figure  C.7:  Pole-zero  plots  for  the  over-modeled  data-adaptive  Toeplitz 
approximation  algorithm. 

Observe  that  as  the  algorithm  converges,  the  poles  (x’s)  are  all  inside  the  unit  circle. 
After  about  50  iterations  their  location  corresponds  to  the  frequency  response  shown 
jr.  Figure  B.14.  Finadly,  we  see  that  the  pole- zero  pairs  different  from  those  shown  in 
Figure  B.ll  are  overlapping  and  tend  to  cancel  each  other,  which  explains  why  both 
of  the  dominant  frequency  components  in  Figure  B.IO  are  present. 
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Figure  C.8;  Over-modeled  example  with  10  dB  of  additive  noise  for  the 
data-adaptive  Toeplitz  approximation  algorithm. 

The  4‘^-order  true  system  was  modeled  with  a  20‘^-order  model.  The  additive  white 
noise  tends  to  slow  down  the  algorithm,  but  only  slightly  distorts  the  frequency 
response  it  is  able  to  achieve.  There  was  10  dB  of  white  noise  added  to  the  output. 
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Figure  C.9:  Pole-zero  plots  for  the  over-modeled  data-adaptive  Toeplitz 
approximation  algorithm  with  lOdB  additive  noise. 

Observe  that  as  the  algorithm  converges,  the  poles  (x’s)  Me  aill  inside  the  unit  circle. 
After  about  100  iterations  their  location  corresponds  to  the  frequency  response  shown 
in  Figure  B.16.  Again,  we  see  that  the  pole-zero  pairs  different  from  those  shown  in 
Figure  B.ll  are  ^dmost  overlapping  and  tend  to  cancel  each  other,  which  explains 
why  both  of  the  dominant  frequency  components  in  Figure  B.IO  are  present. 
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Figure  C.IO:  Parameter  tracks  and  frequency  response  of  the  HR  RLS 
algorithm. 

Observe  that  the  algorithm  converges  to  the  true  solution  and  frequency  response  in 
ten  iterations.  In  this  example  no  noise  was  added  to  the  output,  and  the  performance 
above  is  for  an  exact  modeling  of  the  4**-order  true  system.  Comparing  this  with  the 
data-adaptive  Toeplitz  approximation  algorithm  results  shown  in  Figure  B.IO  reveals 
that  they  behave  differently. 
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Figure  C.ll;  Pole-zero  plots  for  the  IIR  RLS  algorithm. 

Observe  that  as  the  dgorithm  converges  to  the  true  solution,  the  poles  (x’s)  move 
inside  the  unit  circle.  After  about  16  iterations  they  stop  moving,  and  their  location 
corresponds  to  the  frequency  response  shown  in  Figure  B.18.  Finally,  we  see  that 
there  is  jdmost  no  movement  as  the  algorithm  iterates  from  ^  =  16  to  50  iterations, 
which  indicates  that  the  algorithm  has  stabilized. 
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Figure  C.12;  Under-modeled  example  of  the  IIR  RLS  algorithm. 

The  4‘^-order  true  system  was  modeled  with  an  3’’'^-order  model.  The  algorithm 
tries  to  preserve  the  frequency  components  of  the  true  system,  however  it  becomes 
unstable  2dter  about  50  iterations.  No  noise  was  added  to  the  output. 
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Figure  C.13:  Pole-zero  plots  for  the  under-modeled  IIR  RLS  algorithm. 
Observe  that  as  the  algorithm  converges,  the  poles  (x’s)  axe  still  outside  the  unit 
circle.  After  about  50  iterations  their  location  reflects  the  fau:t  that  the  algorithm  is 
unstable,  and  it  fails  to  converge. 
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Figure  C.14:  Over-modeled  example  of  the  HR  RLS  algorithm. 

The  4‘^-order  true  system  was  modeled  with  a  10‘^-order  model.  The  algorithm  also 
fails  to  converge,  as  C2ui  be  seen  after  about  29  iterations  in  the  parameter  track  plot. 
It  is  worth  noting  that  for  a  brief  number  of  iterations,  from  about  A:  =  10  to  27,  the 
algorithm  does  produce  the  true  system  frequency  response.  Unfortunately,  this  is 
not  a  stable  result  of  the  algorithm.  No  noise  was  added  to  the  output. 
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APPENDIX  D:  MATLAB  CODE  FOR  HR 
ADAPTIVE  ALGORITHM 


*/,  This  is  the  HR  ADAPTIVE  algorithm  which  uses  the  Toeplitx  Approximation 
'/,  splitting  R=T+S. 

N  =  input ('How  many  poles  in  the  system  model  would  you  like,  N  =  ?  '); 

M  =  input ('How  many  zeros  in  the  system  model  would  you  like,  M  =  ?  '); 

k  =  input('Enter  the  number  of  iterations  (i.e.  k  =  100).  '); 

noise  =  input ('How  much  noise  would  you  like  (i.e.,  noise  =  .1)  ?  '); 

maxMN  =>  max(M+l,N); 

MN  =  M  +  N; 

s  =  .0000001;  V,  A  small  amount  used  to  start  the  algorithm  below. 

L  =  1; 

P  =  1;  X  Sart  with  one  data  point. 

X  Now  generate  the  "true"  system  data. 

rand ( 'normal') ; 
rand(’seed' ,0) ; 

a  ■  [1,-2.76,3.809,-2.654, .924] ;  X  These  are  the  true  system 
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b  =  .81, -.65, .36] ; 


*L  coefficients. 


X  =  [zeros (M+1 , 1)  ,T3Lnd(P  +  k,l)];  %  The  input  data,  white  noise, 

y  *  f ilter(b,a,x) ;  %  The  output  data. 

X  =  flipud(x) ;  y.  Reorder  the  input  for  convenience, 

rand (' seed' ,5) ;  Y,  A  separate  noise  generator, 

y  =  y  +  sqrt(noise)  *  [zeros  (M+1 , 1)  ;rand(P+k,  1)]  ;  Y,  Add  noise  to  output, 
y  =  flipud(y) ; 

y  ®  [y(l  :P+k,  1)  ;zeros(maxMN+l ,  1)]  ;  Yt  The  additional  zeros  are  used  in  the 
X  ®  [x(l  :P+k,  1)  ;zeros(maxMN+l,  1)]  ;  Y,  computation  of  Z  below,  read  it. 

*/.  Now  I  will  construct  data  matrix,  Z 

Z  *  zeros(P,MN  +  1); 


for  i  =  1:P 

Z(i,:)  =  [y(i+k+2:N+k+l+i)',x(i+k+l:i+k+M+l)']; 

%  The  counting  index 

Y*  ensures  that  the  algorithm  starts  with  zero  data  points, 
end 

%  Now  generate  the  data  correlation  matrix,  R, 

%  and  partition  it  into  4  parts  given  by  M  and  N. 

R  »  Z'  *  Z; 

Ryy  -  R(1:N.1:N); 
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Ryx  =  R(1:N,N+1:M+N+1) ; 

Rxy  =  Ryx ' ; 

Rxx  =  R(N+1:M+N+1.N+1:M+N+1) ; 

%  Now  generate  the  Toeplitz  approximation  to  Ryy  and  Rxx. 

tyy  =  zeros (N, 1) ; 
txx  =  zeros(M+l,l) ; 


for  i  =  1:N 

tyy(i)  =  mean(diag(Ryy,i-l)) ; 
end 

for  i  =  1:M+1 

txx(i)  *  mean(diag(Rxx,i-l))  ; 
end 

Tyy  «  toeplitz(tyy)  +  s  *  eye(N);  %  s  is  added  so  that 

Txx  *  toeplitz(txx)  +  s  *  eye(M+l);  %  these  may  be  inverted. 

Tiy  »  invCTyy) ; 

Tix  *  inv(Txx); 

%  Now  generate  the  cross-correlation  vector,  r, 

X  emd  partition  it  into  two  parts  given  by  M  eind  N. 
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r  =  Z'  *  y(2+k:P+k+l) ; 

ry  *  r(l:N) ; 

rx  =  r(N+l:M+N+l) ; 

*/.  Now  compute  the  initial  estimates  of  the  parameter  vectors. 

b  =  Tiy  *  ry; 
a  =  Tix  *  rx; 

Qk  =  zeros (k,M+N+l)  ;  */,  Holds  the  results  of  each  iteration. 

for  i  =  l:k;  */.  This  is  the  algorithm!!! 

Xn  =  yCk+3-i  ;  k+N+2-i) ;  7,  get  new  output  data 

Xm  *  x(k+2-i  :  k+2+M-i) ;  7,  get  new  input  data 

Yn  *  y(k+2-i) ;  7.  get  next  output  data  point 

%  Use  the  matrix  inverse  lemma  to  compute  the  next 
7.  values  for  Tiy  and  Tix,  note  that  L=l. 

Tiy  “  (Tiy  -  (  Tiy  ♦  Xn  ♦  Xn'  ♦  Tiy)  /  (L+Xn'*Tiy*Xn))/L; 

Tix  *  (Tix  -  (  Tix  *  Xm  ♦  Xm’  ♦  Tix)  /  (L+Xm'*Tix*Xm))/L; 

X  Update  all  of  the  rest  of  the  time-vaxying  quantities. 

Ryy  »  Rjry  +  Xn  *  Xn  ’  ; 
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Ryx  =  Ryx  +  Xn  *  Xm’ ; 

Rxy  ®  Rxy  +  Xm  *  Xn' ; 

Rxx  =  Rxx  +  Xm  *  Xm'; 
ry  =  ry  +  Xn  *  Yn; 
rx  *  rx  +  Xm  *  Yn; 

*/,  Compute  the  next  value  of  the  parameter  vectors. 

b  =  b  +  Tiy  *  (ry  -  Ryx  *  a  -  Ryy  *  b) ; 
a  =  a  +  Tix  *  (rx  -  Rxy  *  b  -  Rxx  ♦  a) ; 

qk(i, :)  =  Cb' ,a'] ; 
end; 
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